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Abstract 

A phase  field  model  is  used  to  numerically  simulate  the  sohdification  of  a pure 
material.  We  employ  it  to  compute  growth  into  an  undercooled  liquid  for  a one- 
dimensional spherically  symmetric  geometry  and  a planar  two-dimensional  rectan- 
gular region.  The  phase  field  model  equations  are  solved  using  finite  difference 
techniques  on  a uniform  mesh.  For  the  growth  of  a sphere,  the  solutions  from  the 
phase  field  equations  for  sufficiently  small  interface  widths  are  in  good  agreement 
with  a numerical  solution  to  the  classical  sharp  interface  model  obtained  using  a 
Green’s  function  approach.  In  two  dimensions,  we  simulate  dendritic  growth  of 
nickel  with  four-fold  anisotropy  and  investigate  the  effect  of  the  level  of  anisotropy 
on  the  growth  of  a dendrite.  The  quantitative  behavior  of  the  phase  field  model  is 
evaluated  for  varying  interface  thickness  and  spatial  and  temporal  resolution.  We 
find  quantitatively  that  the  results  depend  on  the  interface  thickness  and  with  the 
simple  numerical  scheme  employed  it  is  not  practical  to  do  computations  with  an 
interface  that  is  sufficiently  thin  for  the  numerical  solution  to  accurately  represent 
a sharp  interface  model.  However,  even  with  a relatively  thick  interface  the  re- 
sults from  the  phase  field  model  show  many  of  the  features  of  dendritic  growth  and 
they  are  in  surprisingly  good  quantitative  agreement  with  the  Ivantsov  solution  and 
microscopic  solvability  theory. 
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1.  Introduction 


Recently  Kobayashi  [1],  [2],  [3]  reported  computations  of  the  unsteady  phase  field  equa- 
tions in  two  and  three  spatial  dimensions,  which  clearly  showed  the  evolution  of  sohd 
dendritic  structures  into  an  undercooled  melt.  Other  reahstic  features  of  dendritic  growth 
were  also  exhibited  in  these  computations,  such  as  tertiary  side  arm  formation,  side  arm 
coarsening  well  away  from  the  dendrite  tip,  and  the  inclusion  of  hquid  droplets  in  the 
sohd  phase.  This,  to  our  knowledge,  is  the  first  time  that  a computation  of  solidification 
has  been  able  to  show  dendritic  growth  with  the  above  features.  Kobayashi ’s  work  is  a 
qualitative  demonstration  of  the  possible  utihty  of  phase  field  models  of  solidification  as  a 
computational  tool  for  modehng  comphcated,  realistic  solid/liquid  interfaces.  However,  he 
did  not  systematically  investigate  such  important  issues  as  the  accuracy  of  his  numerical 
solutions,  their  relation  to  the  various  classical  formulations  of  solidification,  such  as  the 
Stefan  problem,  or  whether  his  simulations  were  conducted  in  a parameter  regime  that 
corresponded  to  reahstic  growth  conditions  of  an  actual  material. 

It  is  the  focus  of  this  paper  to  address  these  issues,  and  critically  assess  phase  field 
models  as  a viable  computational  technique.  We  go  on  to  compare  the  results  of  our  com- 
putations of  dendritic  growth  using  a phase  field  model  with  current  theories  of  dendrite 
tip  selection,  in  particular,  the  Ivantsov  solution,  marginal  stability  theory,  and  micro- 
scopic solvabihty  theory.  We  also  address  the  issue  of  side  arm  formation.  In  contrast  to 
Kobayashi,  our  work  represents  a quantitative  evaluation  of  phase  field  models,  as  weU  as 
the  theories  of  dendrite  tip  selection  mentioned  above. 

Phase  field  models  of  solidification,  since  their  invention  by  Langer,  see  [4],  following 
an  adaptation  of  the  Model  C proposed  by  Halperin  et  al.  [5],  and  also  independently 
by  Collins  and  Levine  [6],  have  been  subject  to  development  and  rigorous  mathematical 
analysis  by  Caginalp  [7],  [8].  Only  recently  has  computer  technology  advanced  sufficiently 
that  it  is  now  possible  to  numerically  integrate  the  unsteady  phase  field  equations,  in  real- 
istic configurations.  An  early  computation  is  due  to  Smith  [9],  and  more  recently  accurate 
calculations  in  one  spatial  dimension  have  been  obtained  by  Caginalp  and  Socolovsky  [10] 
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and  Lowen  et  al.  [11].  The  advantage  of  the  phase  field  formulation  of  soHdification, 
which  is  described  in  section  2,  is  that  no  distinction  is  made  between  the  sohd,  hquid 
and  interface.  This  allows  the  whole  domain  to  be  treated  in  the  same  way  numerically; 
the  interface  is  not  tracked  but  given  implicitly  by  the  level  set  of  a scalar  function  of 
time  and  space,  the  so-called  phase  field.  In  essence,  the  classical  formulation  of  a free 
boundary  problem  is  replaced  by  a pair  of  nonlinear  reaction-diffusion  equations  for  the 
temperature  and  phase  field.  This  approach  allows  the  computation  of  reahstic  compli- 
cated interfacial  structures  whose  connectedness  changes  in  time.  Numerical  treatments 
of  free  boundary  problems,  using  boundary  integral  or  domain  transformation  methods 
encounter  great  difficulties  in  these  situations. 

Experiments  on  the  free  growth  of  dendrites  by,  for  example,  Glicksman  et  al.  [12], 
show  that  the  dendrite  tip  selects  an  operating  state,  as  characterized  by  the  tip  velocity, 
u,  and  radius  of  curvature,  r,  dependent  on  the  undercooling  of  the  melt,  AT.  Ivantsov’s 
theory  [13]  provides  a local  model  of  a dendrite  tip,  which  is  represented  by  a parabola, 
solidifying  with  constant  velocity  into  an  infinite  melt.  This  theory  determines  the  Peclet 
number,  V = 2ur//c,  as  a function  of  the  undercoohng  parameter,  A = cATIL,  where  c 
is  the  specific  heat,  /c  is  the  thermal  diffusivity,  and  L is  the  latent  heat  per  unit  volume. 
Another  independent  relation  between  v and  r is  required  to  determine  them  uniquely. 
A resolution  of  this  indeterminacy  was  sought  by  introducing  surface  energy  into  the 
theory,  which  introduces  an  additional  length  scale,  do  = cr  j L.,  the  capillary  length  into  the 
problem,  where  a is  the  surface  energy.  This  approach  resulted  in  marginal  stability  theory, 
originally  due  to  Oldfield  [14],  which  predicts  the  other  relation  to  be  C7*  (=  2/cdo/ur^)  = 
0.0192.  This  theory  is  ad  hoc  and  based  on  the  notion  that  the  tip  exists  in  a state  of 
marginal  stability  related  to  the  critical  radius  for  growth  of  a sohd  sphere  or  cylinder  into 
an  undercooled  melt. 

In  contrast,  an  alternate  mathematical  treatment  incorporating  the  surface  energy  is 
known  as  microscopic  solvabihty  theory.  It  is  reviewed  by  Kessler  et  al.  [15]  and  Pomeau 
and  Ben  Amar  [16].  This  theory  indicates  that  surface  energy  in  itself  is  not  sufficient, 
but  that  anisotropy  of  surface  energy  is  required  to  achieve  a steady  stable  tip,  and  cr*  = 
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as  7 — > 0,  where  7 is  the  surface  tension  anisotropy  parameter,  defined  in  two- 
dimensions  by  cr  oc  H-7Cos(fc^),  where  0 is  the  angle  of  the  interface  to  a given  direction  and 
k is  an  integer  that  characterizes  the  symmetry  of  the  crystal.  Unfortunately,  the  physical 
basis  for  this  is  unclear,  the  theory  relying  on  sophisticated  formal  asymptotic  methods.  To 
date,  to  our  knowledge,  there  is  scant  independent  confirmation  of  microscopic  solvabihty 
theory  against  numerical  simulation  or  experiments.  Meiron  [17]  showed  from  a numerical 
boundary  integral  formulation,  that  anisotropy  is  necessary  for  stable,  steady  growth  of 
the  tip  of  a needle  crystal.  Saito  et  al.  [18]  conducted  a numerical  simulation  on  a quasi- 
steady model  using  a Greens  function  technique,  and  were  able  to  roughly  confirm  the 
predictions  of  the  theory.  Experimental  confirmation,  the  ultimate  test  of  such  a theory, 
is  much  less  conclusive.  This  is  because  the  anisotropy  in  the  surface  energy  of  most 
materials  is  difficult  to  measure  and  thus  unknown.  Also  microscopic  solvability  theory 
is  most  weU-developed  for  two  dimensions,  but  experimentally  dendrites  are  usually  three 
dimensional.  Ben  Amar  [19]  has  shown  that  by  choosing  the  best  value  of  the  anisotropy 
parameter  the  experimental  results  of  WilLnecker  et  al.  [20]  could  be  adequately  described 
over  a wide  range  of  values  of  the  dimensionless  undercoofing  parameter.  However,  careful 
experiments  on  pivahc  acid,  camphene  and  succinonitrile  by  Rubinstein  and  Glicksman 
[21],  [22]  do  not  support  microscopic  solvabihty  theory,  in  particular,  they  do  not  confirm 
the  dependence  a*  = as  7 — > 0. 

In  section  2 we  briefly  introduce  a new  phase  field  model  for  sofidification  based  on  an 
entropy  functional  formulation.  In  section  3 we  describe  the  results  of  our  numerical  inte- 
gration of  the  phase  field  equations  in  a spherically  symmetric  geometry,  and  quantitatively 
compare  our  results  to  an  accurate  numerical  solution  of  the  corresponding  free-boundary 
problem  based  on  a Greens  function  technique  due  to  Schaefer  and  GHcksman  [23].  In 
section  4,  we  describe  the  numerical  method  used  to  solve  the  unsteady  phase  field  equa- 
tions in  a planar  two  dimensional  geometry  and  make  a quantitative  assessment  of  it.  We 
compare  our  results  to  both  microscopic  solvability  theory  and  marginal  stability  theory, 
as  weU  as  discuss  side-arm  formation. 

Our  results  in  the  case  of  the  sphericaUy  symmetric  geometry  indicate  that  the  nu- 
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merical  solution  of  the  phase  field  equations  converge  to  the  results  of  the  corresponding 
free-boundary  problem  as  the  interface  thickness  is  reduced  providing  that  the  interface  is 
adequately  resolved.  In  two-dimensional  planar  geometry,  the  validation  of  the  numerical 
solution  of  the  phase  field  equations  is  much  more  difficult  to  perform.  This  is  because 
adequate  resolution  of  the  interfacial  layer  in  a realistic  parameter  range  is  difficult  to 
achieve.  Quahtatively,  our  numerical  solutions  show  dendrite  formation  for  growth  into 
an  undercooled  melt,  and  show  the  importance  of  the  level  of  surface  tension  anisotropy. 
We  find  reasonable  quantitative  agreement  to  the  Ivantsov  solution  and  to  microscopic 
solvabihty  theory  for  a fixed  value  of  the  interface  thickness;  however,  we  show  that  the 
results  are  dependent  on  the  interface  thickness.  This,  we  believe,  provides  the  main  obsta- 
cle to  the  use  of  phase  field  methods  for  the  accurate  computation  of  complex  solid/hquid 
interfaces  for  realistic  growth  configurations  in  more  than  one  space  dimension.  Never- 
theless, the  phase  field  method  produces  many  of  the  quafitative  features  observed  in  real 
crystal  growth,  and  with  further  development  of  the  numerical  solution  techniques  may 
provide  the  best  approach  for  simulating  dendritic  growth. 


2.  The  Governing  Equations 


Phase  field  models  of  the  solidification  of  a pure  material  are  based  on  a Landau-Ginzberg 
free-energy  functional: 


-L 


mT)+-{v^Y 


d^l, 


(1) 


where  Q is  the  region  occupied  by  the  system,  (/>(x,t)  is  the  phase-field,  T(x,  t)  is  the 
temperature  and  e is  a parameter  which  is  constant  for  an  isotropic  material.  The  free- 
energy  density  f[(j>^T)  is  a double  well  with  respect  to  (j).  Various  choices  for  the  precise 
form  of  / have  been  suggested,  the  most  studied  of  which  is 


(2) 


where  a is  a positive  constant,  and  Tm  is  the  melting  temperature  of  the  material.  The 
sohd  and  liquid  phases  are  represented  by  (f)  in  the  neighborhood  of  —1  and  +1  respectively. 
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An  alternative  choice  for  f has  been  proposed  by  Kobayashi  [2]: 

fi<P,T)  = W \-0[T-TM])dC,  (3) 

Jo 

where  14^  is  a constant  and  /5  is  a monotonic  function  of  T — Tm,  with  \/3\  < 1/2  and 
/5(0)  = 0.  This  choice  has  the  advantage  over  Eq.  (2)  that  the  liquid  and  solid  states  are 
represented  by  exactly  equal  to  zero  and  unity  respectively. 

In  both  cases  the  free  energy  Eq.  (1)  is  used  to  derive  a kinetic  equation  for  the  phase 
field  by  requiring  that  it  evolves  in  a manner  such  that  total  free  energy  T decreases 
monotonicaUy  in  time.  The  simplest  choice  is  made  consistent  with  this  requirement: 


(4) 


To  this  is  appended  the  heat  equation  modified  to  take  account  of  the  liberation  of  latent 
heat,  by  the  inclusion  of  an  appropriate  source  term: 


dT 


(5) 


which  provides  an  equation  for  the  temperature,  where  A"  is  a constant  proportional  to  the 
latent  heat  per  unit  volume.  Typically,  Eq.  (5)  is  not  derived  from  basic  thermodynamic 
principles  with  specific  consideration  of  the  form  of  the  free-energy  functional  Eq.  (1), 
so  it  is  not  clear  that  the  solution  of  the  equations  Eq.  (4)  and  Eq.  (5)  will  ensure  that 
J-  is  monotonic  decreasing  in  time.  Penrose  and  Fife  [24]  have  addressed  this  question 
by  employing  the  appropriate  thermodynamic  potential  to  this  non-isothermal  situation, 
namely  an  entropy  functional.  From  this  they  derived,  in  a consistent  manner,  the  phase 
field  equations  when  the  free  energy  density  is  given  by  Eq.  (2).  To  our  knowledge,  the 
phase  field  model  based  on  the  free  energy  density  given  by  Eq.  (3)  has  not  been  placed 
in  a consistent  thermodynamic  setting  in  the  same  way.  However,  this  free  energy  has  the 
advantage  that,  because  the  two  states,  sohd  and  Hquid,  are  given  by  fixed  values  of  (/>, 
0 and  1,  the  latent  heat  released  through  the  source  term  in  the  modified  heat  Eq.  (5)  is 
correctly  accounted  for  in  a numerical  computation  of  the  phase  field  equations  given  by 
Eq.  (4)  and  Eq.  (5).  This  is  not  true  with  the  choice  for  / given  by  Eq.  (2)  as  the  values 
of  (/)  representing  the  sohd  and  hquid  states  depend  upon  e,  T and  a. 
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Another  phase  field  model  which  combines  the  advantages  of  both  the  above  models, 
that  is,  one  that  is  thermodynamically  consistent  and  represents  the  sohd  by  ^ = 0 and 
the  hquid  by  = 1,  has  been  suggested  and  is  discussed  at  length  in  [25 1.  This  new 
model  results  in  the  following  dimensionless  governing  equations  for  the  phase  field  and 
temperature: 


0 — - + 30eQ:Au^(l  — (j)) 


, 1 ^2 


(6) 

(7) 


where  p(^)  = (/5>^(10  — 15<^  + and  prime  denotes  differentiation.  The  solution  of  these 
model  equations  results  in  the  total  dimensionless  entropy  of  the  system  <S,  given  by 


/ 

Jci 


44’, y-)  - 


(8) 


increasing  monotonically  in  time,  with  the  dimensionless  entropy  density  given  by 

44,u)=  f [<(1-C)(C-  l)  + iaAup'{Q]dC  (9) 

Jo 

Here,  length  has  been  scaled  on  some  reference  length  scale  w of  say  the  dimensions 
of  the  domain,  time  on  the  corresponding  thermal  diffusion  time  where  k,  is  the 

thermal  diffusivity,  and  temperature  by  putting  T = Tm  + ATu,  where  AT  is  a reference 
temperature  difference  (such  as,  between  the  melting  temperature  and  the  temperature  at 
the  boundary  of  the  domain). 

The  present  phase  field  model  is  characterized  by  four  dimensionless  parameters.  The 
parameter  A is  the  dimensionless  undercooling,  defined  as 

cAT 


A = 


L 


(10) 


where  c is  the  specific  heat  and  L is  the  latent  heat  per  unit  volume.  In  [25],  it  is  shown  how 
the  remaining  three  constants  are  related  to  the  physical  parameters  which  characterize 
the  interface  dynamics  (i.e.,  interfacial  energy,  cj,  and  mobility,  /z)  and  to  an  estimate  of 
the  interface  thickness,  <f,  which  is  a consequence  of  the  phase  field  approach;  the  following 
definitions  relate  the  remaining  model  parameters  to  the  physical  constants: 

V2wP 


a = 


12caTM  ’ 


(11) 
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(12) 


m = 


kL 


and 


w 


(13) 


Once  the  characteristic  length  scale  w has  been  chosen,  knowledge  of  the  physical  prop- 
erties leaves  one  degree  of  freedom,  namely  e,  which  then  is  used  to  set  the  interface 
thickness.  It  is  expected  that  in  order  to  model  the  physical  behavior  correctly,  the  in- 
terface thickness  must  be  sufficiently  small  compared  to  the  interfacial  macrostructures 
that  we  wish  to  model;  however,  from  a computational  viewpoint,  it  is  desirable  for  the 
interface  thickness  to  be  as  large  as  possible  in  order  that  accurate  solution  of  the  phase 
field  equations  can  be  obtained  for  practical  computational  effort.  This  is  one  of  the  key 
issues  to  be  addressed  by  the  computations  presented  here. 

From  Eq.  (13)  it  is  clear  that  the  hmit  e — ^ 0 corresponds  to  the  interfacial  thickness 
tending  to  zero.  Specifically,  in  [25]  we  show  that  if  we  take  the  hmit  e — ^ 0,  where 
the  constants  a and  m defined  above  are  assumed  to  be  order  one,  we  recover  the  sharp 
interface  model  to  leading  order 


du 

dr 


= 


(14) 


with  interfacial  boundary  conditions 


du 

dn 


liquid 


solid 


a"" 


V, 


nj 


u = -r(— + /C), 


m 


(15) 


(16) 


where  T = crTm/{wLAT)  is  the  dimensionless  capillary  length,  Vn  is  the  dimensionless 
normal  velocity  of  the  interface  (into  the  hquid)  and  JC  is  the  curvature.  The  conservation 
of  heat  across  the  soHd-liquid  interface  allowing  for  latent  heat  production  is  represented 
by  Eq.  (15)  and  the  Gibbs-Thomson  equation  modified  to  account  for  interface  kinetics 
by  Eq.  (16). 
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3.  The  Growth  Of  A Sphere 


Our  first  test  of  the  phase  field  model  is  to  compare  it  to  an  accurate  numerical  solution 
of  the  growth  of  a sphere  into  an  undercooled  melt  originally  described  in  Schaefer  and 
Glicksman  [23].  The  growth  of  the  sphere  is  strongly  influenced  by  curvature,  kinetic,  and 
heat  flow  effects,  thus  providing  the  ingredients  for  a severe  test  of  the  correspondence 
between  sharp  interface  and  phase  field  descriptions  in  the  limit  e — > 0. 

Schaefer  and  Glicksman  [23]  solved  the  sharp  interface  problem  Eq.  (14),  Eq.  (15), 
and  Eq.  (16),  recovered  in  the  Hmit  e — > 0 in  a spherically  symmetric  geometry  using 
a Greens  function  method.  We  have  used  the  same  technique  here  modifying  slightly 
their  normahzation  to  be  consistent  with  that  we  use  here  in  the  phase  field  model.  The 
initial  data  corresponded  to  a solid  sphere  of  radius  po  chosen  to  be  10%  greater  than  the 
equilibrium  radius  p* . In  our  non-dimensionalization,  we  choose  p*  as  the  reference  length 
scale  ly,  and  AT  to  be  equal  to  the  undercooHng  (taken  to  be  positive).  Thus,  initially 
we  set  u = —1  everywhere.  The  dimensionless  equihbrium  radius  is  given  from  Eq.  (16) 
as  p*  = 1/(3\/2q;A),  and  because  we  choose  p*  — \ this  gives  a = l/(3\/2  A).  Also,  we 
introduce  the  parameter 


PlcfTm 
2k,L  ’ 

which  represents  the  mobihty  of  the  interface  and  is  the  same  dimensionless  group  em- 
ployed in  [23],  and  we  note  that  m = 2^.  The  nondimensional  boundary  condition  for  the 
interfacial  temperature  Eq.  (16)  is  given  in  terms  of  ^ as 


dp 

dr 


-4^ 


(18) 


3.1.  Numerical  Method  For  Solution  Of  the  Sharp  Interface  Problem 


As  the  sphere  grows,  heat  is  emitted  from  its  surface  at  a rate  proportional  to  the  latent 
heat  and  to  the  growth  rate.  The  resulting  temperature  field,  and  in  particular  the  surface 
temperature  of  the  growing  sphere  at  r = /?('r),  is  calculated  by  a Green’s  function  integral: 

= + ^ G{p,p'{r'),T,T')^dT',  (19) 
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where  for  a sphere 


G{ 


r,  r , T,  T 


') 


( jr  - r'Y  \ 

\ 16(t  — r'Y  ) 


— exp 


The  radius  of  the  sphere  is  given  simply  by 


16(t  — t'Y  ) 


(20) 


The  numerical  procedure  used  to  increment  the  interface  position  pn  and  temperature 
at  the  nth  time  step  {r  = Tn  = nAr)  to  the  n+lth  time  step  (r^+i  = + At)  was  as 

follows: 


1.  Compute  dp! dr  at  r = r^,  denoted  by  pn,  from  Eq.  (18). 

2.  Compute  p^+i  from  Eq.  (21)  employing  a backward  difference,  i.e.,  pn+i  = pn-^PnAr. 

3.  Update  the  interfacial  temperature  from  Eq.  (19)  by 

i=n 

u{pr,,nAr)  = At  p^.iG{pr„  p^,nArYAr). 

t=i 

An  alternative  quadrature  formula  of  the  Greens  function  integral  in  step  3 was  also  tried 
using  only  pi_i  on  each  time  interval  and  was  found  to  converge  much  less  rapidly.  For  the 
conditions  considered  in  this  paper,  calculations  were  carried  out  using  a range  of  values 
of  the  time  step.  It  was  found  that  for  time  steps  smaller  than  a value,  which  depends  on 
A and  the  results  agreed  within  0.1%. 


3.2.  Numerical  Method  For  The  Phase  Field  Model 

We  assume  that  initially  the  temperature  is  everywhere  uniform  at  u = —1.  We  employed 
a finite  difference  scheme  to  solve  the  phase  field  equations  Eq.  (6)  and  Eq.  (7)  expressed 
in  a spherically  symmetric  geometry  in  which  the  sole  spatial  variable  is  r,  the  radial 
distance.  Second-order  central  differences  were  used  to  discretize  the  Laplacian  operator 
on  a uniform  spatial  mesh  for  the  domain  0 < r < R,  with  AT-f  1 nodes.  Unhke  the  Greens 
function  approach  described  above  the  phase  field  computation  requires  a finite  domain 
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and  thus  we  must  be  careful  in  choosing  appropriate  outer  boundary  conditions  a-t  r = R. 
We  chose  them  to  be 


du 

dr 


li  + 1 

R 


= 0, 


and 


d(l) 

dr 


= 0 at  r = J?, 


(22) 


which  is  consistent  with  assuming  a quasi-static  approximation  for  the  heat  diffusion  equa- 
tion in  the  far  field.  The  phase  field  is  simply  equal  to  unity  in  the  far  field.  At  the  origin 
the  coordinate  singularity  in  the  Laplacian  operator  was  avoided  by  employing  a local 
form  for  the  solution  as  described  in  Smith  [9]  and  Neumann  boundary  conditions.  The 
solution  was  advanced  in  time  using  a second-order  Crank-Nicholson  scheme  for  the  Lapla- 
cian terms  and  a backward  difference  for  the  latent  heat  production  term  in  the  modified 
heat  equation.  Specifically,  the  solution  was  advanced  a single  time  step  in  the  following 
manner: 

1.  Update  the  phase  field  by  solving  the  tridiagonal  system, 


2.  Update  the  temperature  by  solving  the  tridiagonal  system. 


Here  and  are  the  discrete  approximations  to  ^ and  u,  etc  represent  the  tridi- 

agonal matrices  that  result  from  the  temporal-spatial  discretization, 
and  !,($-)  - The  computations  were  conducted  on  a Cray-YMP  and  the 

solution  of  the  tridiagonal  systems  was  conducted  using  the  Cray  library  routine  TRID 
which  employs  cyclic  reduction.  The  code  was  well  vectorized  and  achieved  in  excess  of 
200  MFLOPS  on  a single  processor. 

3.3.  Comparison  of  the  Two  Methods 

We  conducted  calculations  using  both  methods  for  an  undercooling  of  0.5  and  two  different 
values  of  the  mobility,  ( = 0.05  and  1.0.  The  most  comprehensive  set  of  computations 
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was  conducted  for  = 1.0.  For  this  case  in  conducting  the  computations  on  the  phase 
field  model,  we  used  five  different  values  of  e (0.1,  0.075,  0.05,  0.02  and  0.01)  and  five 
different  grids  with  mesh  spacings  Ar  (0.001,0.002,0.01,0.02  and  0.1)  for  each  value  of  e. 
The  time  step  was  taken  to  be  lO""^.  It  was  found  that  there  was  a critical  value  of  the 
time  step,  that  was  apparently  insensitive  to  the  mesh  size,  above  which  the  numerical 
scheme  was  unstable,  despite  the  impficit  discretization  of  the  diffusion  operator.  We 
ascribe  this  to  the  explicit  discretization  of  the  source  term^  in  the  phase  field  and  heat 
equation.  It  was  found  that  invoking  the  outer  boundary  condition  at  r = 10  was  sufficient 
for  approximating  an  unbounded  domain. 

For  the  case  ( = 0.05,  which  is  a more  reahstic  value  for  a supercooled  pure  metal, 
calculations  were  only  performed  on  the  phase  field  model  for  e = 0.01,  with  Ar  = 0.05, 
and  At  = 10“^,  for  which  we  were  able  to  obtain  excellent  agreement  with  the  numerically 
computed  solution  of  the  classical  formulation,  but  with  the  outer  boundary  placed  further 
away  at  r = 20.  We  attribute  this  to  the  increased  communication  by  the  temperature 
field  to  the  far-field  in  this  case  of  a lower  interfacial  mobihty  w_.  'e  the  interface  motion 
is  relatively  sluggish  relative  to  the  thermal  diffusion  time.  Henceforth,  we  restrict  our 
discussion  to  the  case  = 1.0. 

On  the  finest  grids  excellent  agreement  was  obtained  between  the  two  procedures  as 
shown  in  Figure  1 for  Ar  = 0.001  and  e = 0.01.  It  was  found  that  the  numerical  solution 
exhibited  a temporal  oscillation  when  the  mesh  was  too  coarse.  In  all  except  one  of  these 
cases  the  breakdown  corresponded  to  values  of  e/Ar  (a  measure  of  the  resolution  of  the 
interface  by  the  computational  mesh)  less  than  or  equal  to  unity.  In  particular,  for  a fixed 
value  of  e,  on  increasing  the  mesh  spacing  the  numerical  instabihty  first  manifested  itself 
as  a significant  oscillation  in  the  interface  velocity;  at  larger  values  of  Ar  it  was  apparent 
as  a significant  oscillation  in  the  surface  temperature  also.  These  observations  broadly 
confirm  the  suggestion  by  Osher  [26]  and  Caginalp  and  Socolovsky  [10]  that  the  numerical 
solution  of  the  phase  field  equations  by  a finite  difference  method  on  a uniform  mesh  could 
be  expected  to  breakdown  for  Ar  > e. 

We  now  discuss  the  solutions  computed  on  sufficiently  fine  meshes  for  which  the  numer- 
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Position,  r 


Figure  1:  Tlie  solid  curves  (from  bottom  to  top)  are  the  temperature  at  times  r = 1.25, 
2.5,  3.75,  . . . , 11.25,  as  a function  of  the  radial  coordinate,  computed  from  the  phase  field 
model  for  the  case  e = 0.01,  A = 0.5  and  ^ = 1.0.  The  dashed  curve  represents  the  locus 
of  the  interface  (j)  = 0.5  and  the  almost  coincident  solid  circles  are  taken  from  the  Greens 
function  solution  at  the  corresponding  time  levels. 
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At 

e =0.01 

e =0.02 

e =0.05 

e = 0.075 

e = 0.1 

0.001 

1.520 

-0.7502 

-3.347 

-5.213 

-7.058 

0.002 

1.504 

-0.7540 

-3.347 

-5.213 

-7.058 

0.01 

1.009 

-0.8753 

-3.361 

-5.2173 

-7.060 

0.02 

92.44 

-1.2452 

-3.404 

-5.2300 

-7.063 

0.1 

91.49 

9.149 

-3.086 

-5.6091 

-7.063 

1 

Table  1:  The  relative  difference  in  the  surface  temperature  calculated  by  the  two  methods 
over  an  interval  of  12  time  units  expressed  as  a percentage.  The  entries  in  italics  indicate 
that  the  solution  was  subject  to  a numerical  oscillatory  instabiHty.  Those  in  bold  were 
employed  in  the  curve  fitting  discussed  in  section  3.3;  in  particular,  the  emboldened  entries 
in  the  columns  were  fitted  as  —0.748337  — 0.00153s  — 0.52s^,  for  e = 0.02;  —3.34638  — 
0.000375s  - 0.36s2,  for  i = 0.05;  5.21304  - 0.00462s  - 0.22s2,  for  i = 0.075;  -7.05878  - 
0.0330495s  + 0.284232s^,  for  e = 0.1,  where  s = Ar/e. 

ical  oscillation  was  not  present.  We  found  the  best  agreement  in  the  interface  temperature 
between  the  two  methods  was  obtained  when  it  was  expressed  as  a function  of  the  in- 
terface position,  the  disparity  being  insensitive  to  e/Ar  (in  the  phase  field  formulation 
the  interface  position  was  defined  to  be  given  by  = 1/2)-  The  interface  temperature 
as  a function  of  interface  position  is  controlled  largely  by  heat  flow,  and  this  excellent 
agreement  indicates  that  our  scheme  provides  a good  solution  of  the  heat  equation.  How- 
ever, the  agreement  between  the  two  methods  for  the  interface  temperature  Ui  (as  well 
as  its  position  p and  velocity  p),  expressed  as  a function  of  time  depended  much  more 
strongly  on  the  values  e and  Ar.  In  particular,  the  disparity  between  the  two  numerical 
methods  increased  monotonically  with  time.  This  distinction  between  the  errors  in  the 
interface  temperature  expressed  as  a function  of  time  or  position  is  illustrated  in  Figure 
2 for  e — 0.1.  The  central  dashed  line  indicates  the  interface  temperature  as  a function 
of  position,  as  computed  by  the  phase  field  method.  This  curve  is  almost  coincident  with 
the  solid  circles  which  represent  the  same  quantity  but  obtained  from  the  Greens  function 
method.  However,  each  solid  circle  is  not  coincident  with  the  interface  temperature  for 
the  corresponding  time  level,  thus  indicating  a considerable  disparity  in  the  two  methods 
when  the  interfacial  temperatures  are  compared  as  a function  of  time. 
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Position,  r 


Figure  2:  The  solid  curves  (from  bottom  to  top)  are  the  temperature  at  times  r = 1.25, 
2.5,  3.75,  . . . , 11.25,  as  a function  of  the  radial  coordinate,  computed  from  the  phase  field 
model  for  the  case  e = 0.1,  A = 0.5  and  { = 1.0.  The  dashed  curves  are  the  locus  of  (from 
left  to  right)  cj)  = 0.1,  (j)  = 0.5  and  (j)  = 0.9.  The  solid  circles  represent  the  locus  of  interface 
temperatures  computed  from  the  Greens  function  solution;  each  sofid  circle  corresponds 
to  the  same  time  levels  as  the  solid  curves,  and  the  time  levels  increase  upwards. 
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A summary  of  the  errors  between  the  two  methods  is  given  Table  1,  where  the  error 
between  the  two  techniques  in  the  interface  temperature  (averaged  as  a function  of  time 
over  the  computational  time  interval  of  12  time  units)  is  tabulated  in  terms  of  the  mesh 
spacing  and  e.  The  itaHcized  entries  indicate  an  unstable  numerical  solution,  as  discussed 
above.  The  variation  in  the  errors  of  p and  p showed  similar  behavior  to  the  interface 
temperature.  For  a fixed  value  of  the  mesh  spacing  the  error  between  the  two  techniques 
decreases  as  e decreases,  except  for  the  case  of  reduction  of  e from  2x10“^  to  the  smallest 
value  10“^.  This  we  attribute  to  the  error  between  the  exact  phase  field  solution  and 
the  classical  solution  being  sufficiently  small  for  e = 0.01  that  it  is  comparable  to  the 
numerical  errors  in  the  Greens  function  solution  of  the  classical  problem,  as  weU  as  the 
errors  in  the  quadrature  formula  (with  a time  step  of  0.05)  used  to  compute  the  average 
surface  temperature  from  the  numerically  computed  phase  field  solution.  Linear  regression 
of  the  errors  against  e for  fixed  values  of  Ar  show  the  error  between  the  numerical  solution 
of  two  different  formulations  to  be  first  order,  i.e.,  oc  Ae,  where  the  constant  A was  found  to 
be  78.6,78.56,77.13,73.32  for  Ar  = 0.001,0.002,0.01,0.02,  respectively.  The  convergence 
with  i is  therefore  insensitive  to  the  mesh  size  (assuming  that  numerical  instabihty  is  not 
present).  It  also  confirms  the  asymptotic  theory  of  the  phase  field  equations  in  the  limit 
e — > 0 by  Caginalp  [8]  which  indicates  that  their  solution  approaches  the  sharp  interface 
solution  in  this  limit.  This  strong  hnear  dependence  of  the  error  on  e for  a fixed  mesh 
allows  the  application  of  accelerated  convergence  techniques,  such  as  Richardsons  method, 
to  the  phase-field  equations. 

With  the  value  of  e fixed  the  errors  decreased  as  the  mesh  spacing  was  refined,  albeit 
very  slowly.  There  was  one  exception  to  this  for  e = 0.01,  where  regression  of  the  errors 
with  respect  to  Ar  indicates  that  this  is  not  so.  This  we  attribute,  as  discussed  above, 
to  comparable  errors  associated  with  the  Greens  function  technique  and  the  quadrature 
formula  used  when  i is  sufficiently  small.  We  fitted  a quadratic  polynomial  to  the  errors 
as  a function  oi  s = Ar/e  for  e = 0.02,  0.05,  0.075  and  1.0.  The  values  of  the  coefficients 
are  given  in  Table  1 and  show  that  the  convergence  is  predominantly  quadratic;  however, 
the  coefficients  do  depend  on  e. 
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In  summary,  these  calculations  reveal  that  a numerical  instability  occurs  when  either 
the  time  step  exceeds  a critical  value,  which  is  independent  of  the  mesh  size,  or  the  mesh 
spacing  is  approximately  greater  than  or  equal  to  e.  The  optimum  value  of  Ar  is  given  by 
Ar  ^ e because  the  error  is  not  significantly  reduced  by  using  a finer  mesh  compared  to 
the  consequent  increase  in  processor  time.  The  errors  themselves  appear  to  be  0{e)  and 
dependent  only  weakly  on  the  mesh  spacing  provided  Ar  < e. 


4.  Two-Dimensional  Calculations 


In  addition  to  the  spherically  symmetric  calculations  previously  described,  a series  of  two- 
dimensional  numerical  calculations  are  performed  in  order  to  further  evaluate  the  present 
phase  field  model.  The  primary  objective  is  to  evaluate  the  behavior  of  the  given  phase  field 
model  for  reahstic  physical  parameter  values,  i.e.,  evaluate  how  well  the  diffuse  interface 
model  reproduces  the  essential  physical  mechanisms,  and  investigate  the  feasibility  of 
accurate  numerical  simulations  of  unconstrained  solidification  of  a pure  material  in  an 
undercooled  melt. 

For  the  two-dimensional  calculations,  the  phase  field  equation  given  by  Eq.  (6)  is 
modified  to  include  anisotropy  in  the  parameter  e,  where  the  anisotropy  is  assumed  to 
have  the  form 

e(^)  = eT){9)  = e(l  + 7 cos  kO). 


The  angle  6 is  defined  as  the  angle  between  the  normal  to  the  interface  and  the  x-axis, 
and  k specifies  the  mode  number.  Evaluating  the  variational  derivative  of  Eq.  (1)  after 
including  the  variation  of  e with  orientation  yields  a modified  version  of  Eq.  (6): 


dcf) 


(f>  — - + 30eaAu(^(l  — <^) 


— e 


-2 


dx 


ivWv'iS) 


dy. 


■ (r?^(e)V^)| . (23) 

In  the  phase  field  model,  the  interface  is  represented  by  level  sets  of  (j).  In  order  to  compute 
the  anisotropic  behavior,  the  orientation  angle  is  determined  in  terms  of  the  phase  field  cj) 
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using  the  following  relation  for  the  normal  vector: 


iv^l 

From  this  expression  we  have  the  definition 


n = = cos  9x  + sin  dy. 


tan^  = — , 


and,  in  addition,  we  obtain 


n 4^x4^xy  4’y^xx  ^ 4^y4^xy 

\V(j)\^  ’ ''''  ■ 

These  relations  are  used  to  compute  the  terms  in  Eq.  (23)  which  arise  from  expanding  the 
derivatives  of  e{6). 

It  is  shown  in  [27]  that  the  only  modification  to  the  free  boundary  problem  Eq.  (14)  to 
Eq.  (16),  obtained  in  the  limit  i ^ 0,  is  to  the  modified  Gibbs-Thomson  equation  which 
becomes 

We  note  that  allowing  e to  depend  on  9 modifies  both  the  interface  kinetic  and  curvature 
terms.  It  is  shown  in  [27]  that  the  latter  is  the  same  as  that  obtained  from  sharp  interface 
models  where  the  surface  energy  depends  on  the  interface  orientation. 

For  the  two-dimensional  simulations,  symmetry  conditions  (vanishing  Neumann  con- 
ditions for  both  (j)  and  u)  are  applied  at  the  boundaries  of  the  rectangular  domain,  which 
has  scaled  dimensions  of  Xl  and  Yl  in  the  x and  y coordinate  directions,  respectively. 
We  note  that  the  results  presented  here  for  the  temperature  and  phase  field  are  displayed 
after  reflecting  the  computational  domain  about  the  hne  y = 0,  which  corresponds  to  the 
axis  of  the  dendrite.  From  the  chosen  definition  of  the  dimensionless  variables,  the  value 
u = 0 corresponds  to  the  melting  temperature  of  the  pure  material,  while  u = —1  is  the 
undercooling  temperature  (the  dimensional  undercooling  is  AT).  Initially  for  each  calcu- 
lation, a small  region  of  solid  (u  = 0)  is  located  at  the  x = 0,  y = 0 corner  of  the  domain, 
and  the  remainder  of  the  domain  is  undercooled  liquid  (u  = ~1)-  The  shape  of  the  initial 
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solid  region  is  either  one-quarter  of  a circle  with  an  initial  radius  denoted  by  or  one- 
quarter  of  an  elhpse  with  semi-minor  axis  i/o  and  semi-major  axis  Xo-  Either  continuous 
initial  conditions,  for  which  the  values  of  (f)  and  u change  smoothly  over  a thin  region,  or 
discontinuous  conditions  were  used  in  the  computations.  For  the  results  presented,  the 
long-time  behavior  of  the  solutions  is  insensitive  to  the  type  of  initial  conditions  assumed. 

The  two-dimensional  simulations  were  performed  using  property  values  for  pure  nickel. 
The  values  for  the  required  material  parameters  are:  a = 3.7x10“^  J/  cm^,  = 1728  K, 
L = 2350  J/  cm^,  c = 5.42  J/  K cm^,  /c  = 0.155  cm^/s,  and  /x  = 285  cm/K  s.  Except 
for  the  kinetic  coefficient,  /i,  these  property  values  for  nickel  are  readily  available  (see 
[28]).  The  value  for  the  kinetic  coefficient  lies  in  the  range  of  estimated  values  and  was 
chosen  here  partly  from  consideration  of  the  numerical  values  of  the  dimensionless  model 
parameters. 

In  order  to  completely  determine  the  dimensionless  parameters  given  by  Eq.  (11)  - 
Eq.  (13),  we  must  choose  values  for  the  reference  length,  u;,  and  the  interface  thickness 
6.  The  choice  of  these  parameters  is  based  on  the  physical  structure  we  wish  to  compute 
and  the  practical  Hmitations  of  accurately  resolving  gradients  within  the  interfacial  region 
for  a desired  computational  domain  of  size  Xl  and  Yl-  For  the  simulation  of  dendritic 
growth,  we  choose  to  relate  8 and  w to  an  estimate  of  the  dendrite  tip  radius.  Clearly, 
we  would  expect  that  8 must  be  much  smaller  than  the  tip  radius  of  the  dendrite  and 
that  the  domain  must  be  many  times  larger  than  the  tip  radius  to  simulate  the  growth 
and  development  of  the  dendrite.  In  order  to  begin  the  calculations,  we  estimated  the 
dendrite  tip  radius  based  on  marginal  stability  theory  [14]  for  a given  value  of  undercoohng, 
represented  nondimensionally  here  by  A.  The  computations  are  facihtated  by  lower  values 
of  A (large  values  of  undercooling,  AT),  because  the  dendrite  grows  faster  and  encompasses 
a larger  portion  of  the  computational  domain.  Values  of  A in  the  range  of  0.4-0. 5 are  used 
here.  A dimensionless  undercoohng  of  A = 0.5  corresponds  to  an  actual  undercoohng  of 
217  K for  nickel,  which  is  an  attainable  level  of  undercoohng  [20].  For  an  undercoohng 
corresponding  to  A = 0.5,  marginal  stabihty  gives  an  estimated  value  for  the  tip  radius 
of  1.7  X 10~^  cm.  Based  on  some  prehminary  computational  experiments  to  determine 
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how  the  interface  thickness  and  the  size  of  the  domain  should  be  chosen  in  relation  to 
the  tip  radius,  we  obtained  the  parameter  value  a — 400  in  the  definition  Eq.  (11).  For 
the  physical  parameters  of  nickel  given  above,  this  yields  a value  for  the  length  scale 
w = 2.1  X 10 cm;  the  dimensionless  parameter  m has  the  value  0.05  for  nickel. 

With  the  parameters  a and  m specified,  the  parameter  e was  used  in  the  computations 
to  vary  the  thickness  of  the  interface.  The  computational  resolution  was  determined  by 
the  extent  of  the  domain  chosen  [Xl  and  Yl)  and  the  number  of  grid  points  used  in  the 
discretization  of  the  domain.  We  performed  computations  with  e equal  to  0.005,  0.0033, 
and  0.0025,  but  for  the  majority  of  the  simulations  of  dendritic  growth  into  an  undercooled 
liquid  we  used  the  value  e = 0.005. 

4.1.  Numerical  Method 

One  of  the  clear  advantages  of  the  phase  field  model  is  that  the  location  of  the  solid/hquid 
interface  does  not  have  to  be  determined  expHcitly.  However,  accurate  numerical  solutions 
to  the  phase  field  equations  require  that  gradients  of  the  field  variables  be  adequately 
resolved  over  the  thin  interfacial  region.  Our  objective  is  to  evaluate  the  behavior  of  the 
phase  field  model  presented  here  by  simulating  the  growth  of  a two-dimensional  dendrite 
into  an  undercooled  liquid.  In  order  to  perform  the  simulation,  the  governing  equations 
are  solved  numerically  in  a two-dimensional  rectangular  domain  of  sufficient  size  to  allow 
for  the  development  of  characteristic  dendritic  structure  (e.g.,  side  arms).  Even  with  the 
phase  field  approach,  an  optimum  solution  procedure  should  employ  some  type  of  adaptive 
technique,  particularly  for  the  phase  variable,  since  the  solution  varies  over  such  a small 
part  of  the  domain.  Since  our  interest  is  primarily  to  evaluate  the  model  itself,  we  have 
chosen  to  use  straightforward  finite  difference  solution  techniques  applied  on  a uniform 
computational  grid.  An  advantage  of  this  approach  is  that  the  numerical  implementation 
of  such  techniques  is  quite  simple  and  it  is  easy  to  take  full  advantage  of  highly  vectorized 
large-scale  computers. 

The  governing  equations  given  by  Eq.  (7)  and  Eq.  (23)  are  a pair  of  coupled,  second- 
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order,  nonlinear  parabolic  equations.  They  are  discretized  spatially  using  second-order 
finite  differences  on  a uniform  grid  characterized  by  mesh  spacings  AX  and  AF  in  the 
X and  y coordinate  directions,  respectively;  for  the  temporal  discretization  we  introduce 
the  time  step  At.  In  order  to  maximize  the  computational  efficiency,  we  employ  explicit 
time-differencing  on  the  (j)  equation  which  is  nonlinear  in  all  terms  except  for  the  highest- 
order  spatial  derivatives.  The  heat  Eq.  (7)  is  linear  in  the  temperature  u but  contains  the 
source  term  depending  on  (j).  With  explicit  time- differencing  Eq.  (7)  would  be  subject  to 
a more  restrictive  time  step  requirement  than  Eq.  (23);  thus,  we  employ  the  alternating- 
direction  implicit  method  (ADI)  on  Eq.  (7)  which  for  the  simple  linear  heat  equation 
is  unconditionally  stable  and  second-order  accurate  in  space  and  time.  The  methods 
employed  here  are  described  in  many  standard  texts  on  finite  difference  techniques  for 
partial  differential  equations  (for  example,  Richtmyer  and  Morton  [29]). 

For  the  nonhnear  phase  field  equation  with  exphcit  time- differencing  (e.g.,  simple  Euler, 
Adams-Bashforth),  the  unsteady  linear  diffusion  part  of  the  equation  is  subject  to  the 
stability  restriction:  At  < (AX)^/(4m);  however,  the  additional  nonhnear  terms  in  the 
equation  may  impose  a more  restrictive  condition  on  the  size  of  the  allowable  time  step. 
The  actual  size  of  the  time  step  was  determined  by  numerical  experimentation.  Both  first- 
and  second-order  accurate  explicit  time  differencing  were  evaluated  for  solving  Eq.  (23); 
however,  the  first-order  method  was  used  to  obtain  the  results  presented  here,  since  the 
restriction  on  the  size  of  the  time  step  for  the  (j)  equation  was  such  that  first-order  accuracy 
in  time  proved  to  be  adequate  to  obtain  a good  balance  in  the  spatial  and  temporal 
truncation  error. 

4.2.  Computation  of  Dendrites 

Using  the  numerical  technique  described  above  we  have  sought  to  compute  the  evolution 
of  dendritic  structures  using  the  parameter  values  m = 0.05  and  a = 400  which  were  deter- 
mined from  the  physical  properties  for  nickel  as  described  previously.  The  dimensionless 
parameters  e.  A,  and  7 are  varied  in  the  calculations  as  is  the  size  of  the  computational 
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domain.  We  assume  4-fold  anisotropy  for  all  the  calculations  presented  here  so  the  mode 
number  k has  the  value  4. 

We  found  that  in  the  absence  of  anisotropy  (7  = 0)  the  solid  developed  as  an  amorphous 
structure  in  which  the  dendrite  tips  were  subject  to  repeated  tip  spHtting.  A typical 
example  of  this  behavior  is  shown  in  Figure  3(a).  For  nonzero  values  of  the  anisotropy 
parameter  7,  a very  distinctive  needle  crystal  formed  along  each  coordinate  direction 
as  a result  of  the  four-fold  symmetry  of  the  anisotropy  as  shown  in  Figure  3(b).  This 
is  consistent  with  the  recent  microscopic  solvability  theory  of  dendrite  growth,  which 
indicates  the  importance  of  surface  energy  anisotropy,  as  well  as  calculations  on  local 
models  of  solidification  [15]. 

For  smaller  values  of  the  anisotropy  (7  less  than  approximately  0.01),  the  tip  did  not 
settle  down  to  a steady  state  over  the  period  of  the  computation;  it  is  not  clear  from 
our  calculations  whether  a steady  state  would  have  occurred  on  a larger  computational 
domain  over  a longer  period  of  dimensionless  time.  However  it  is  clear  that  anisotropy  has 
a profound  effect  on  the  evolution  of  the  crystal.  For  values  of  the  anisotropy  parameter 
10“^  < 7 < 2 X 10~^  the  dendrite  tip  rapidly  locked  into  a definite  steady  operating  state 
for  the  computational  domains  we  employed.  Unless  otherwise  stated  aU  computations 
presented  below  for  quantitative  comparison  with  non-zero  anisotropy  yielded  a needle 
crystal. 

In  order  to  quantitatively  describe  dendritic  growth  we  compute  the  the  tip  temper- 
ature U,  radius  5,  and  velocity  V ( overbars  denote  phase  field  results)  for  the  dendrite 
aligned  with  the  x-axis  in  the  following  manner:  we  define  the  interface  by  the  locus 
(/)(a:,y,t)  =1/2  and  so  estimate  the  the  dendrite  tip  position  on  the  x-axis,  and  the  corre- 
sponding temperature  and  curvature,  at  each  time  level,  by  linear  interpolation  from  the 
mesh.  The  radius  of  curvature  of  the  tip  was  approximated  at  each  mesh  point  on  the 
x-axis  by  employing  the  identity 


^ ^4>yy 
R 4>x  ’ 

where  second-order  accurate  finite  differences  were  used  to  approximate  the  partial  deriva- 
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Figure  3:  The  computed  phase  field  for  parameters  corresponding  to  nickel  with  e = 0.005, 
A = 0.5,  Xl  = 4.5  and  Yl  = 2.25;  note  that  the  computational  domain  has  been  reflected 
about  the  horizontal  centerhne.  The  results  show  the  the  effect  of  anisotropy  level: 

(a)  7 = 0 and  (b)  7 = 0.01.  The  black  area  is  solid  0 < ^ < 0.1;  the  white  area  is  the 
interfacial  region  O.l  < (f)  < 0.9;  the  gray  area  is  hquid  0.9  < (f)  < 1. 
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tives.  The  tip  velocity  was  approximated  at  the  the  midpoint  of  each  time  interval  by  using 
central  differences  on  the  dendrite  tip  position. 

In  Figure  4 we  plot  V,  R and  as  a function  of  time.  There  is  clearly  a very  well- 
defined  operating  state  for  the  tip  over  the  greater  part  of  the  time  interval,  following  an 
initial  transient  as  the  tip  quickly  attains  its  operating  state.  The  final  transient  is  due  to 
the  dendrite  tip  encountering  the  end  of  the  box. 

In  performing  the  numerical  computations  several  issues  had  to  be  addressed  in  order 
to  evaluate  the  behavior  of  the  present  phase  field  model.  First,  we  consider  the  effect 
of  spatial  mesh  size  and  time  step  for  a fixed  set  of  dimensionless  parameters.  In  Table 
2,  we  list  computed  values  of  the  interface  velocity  and  dendrite  tip  radius  for  different 
mesh  spacings  and  time  steps.  We  note  here  that  the  mesh  spacing  in  each  direction,  AX 
and  Ay,  are  always  taken  to  be  equal  to  one  another  in  all  the  calculations.  For  the  finer 
spatial  meshes,  the  time  step  restriction  of  the  explicit  treatment  of  the  (f)  equation  only 
permitted  stable  numerical  solutions  for  the  smaller  time  step  values.  It  is  apparent  from 
the  results  that  the  tip  velocity  is  more  sensitive  to  the  discretization  error  than  the  tip 
radius.  We  note  that  the  computation  time  for  AX  = 0.0025  and  At  = 2.5  x 10“^  is 
16  times  greater  than  for  the  case  with  AX  = 0.005  and  At  = 1 x 10“^,  but  with  only 
a corresponding  change  of  about  10  % in  the  numerical  values.  In  fight  of  this,  for  the 
calculation  of  dendrites  presented  below,  we  employ  the  coarser  of  these  two  meshes  and 
the  larger  of  the  time  steps  in  order  to  adequately  investigate  the  parameter  space.  In  the 
calculations  presented  here,  we  estimate  that  there  are  approximately  nine  spatial  mesh 
points  within  the  interfacial  region  when  AX  = e;  this  level  of  resolution  was  used  unless 
otherwise  stated.  We  note  that  in  a similar  computation,  Kobayashi  [2]  employed  a mesh 
approximately  six  times  coarser. 

We  now  discuss  the  effect  of  the  interface  thickness  on  the  computed  results.  For  the 
parameter  values  corresponding  to  nickel  (i.e.,  a = 400  and  m = 0.05),  an  interface  thick- 
ness (defined  to  be  represented  by  (j)  values  between  0.1  and  0.9)  corresponding  to  e = 
0.005  was  found  to  be  the  largest  allowable;  for  e equal  to  0.01,  the  interfacial  region  is 
spread  apart  to  such  an  extent  that  it  no  longer  approximates  a thin  interface,  while  for 
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Figure  4:  The  estimated  values  of  the  dimensionless  dendrite  tip  velocity  V,  temperature 
U,  and  radius  R,  against  time  r for  the  case  e = 0.005,  7 = 0.015,  and  A = 0.4. 
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AX  = 0.003 

V R(IO^) 

AX  = 0.0025 

V H(102) 

1.0  xlO-"* 

5.0  xlO-^ 
2.5  xlO-^ 

4.57  7.34 

4.70  7.14 
4.84  7.11 

5.43  6.79 

5.62  6.76 

5.72  6.74 

5.85  6.65 
5.96  6.61 

6.04  6.53 

Table  2;  The  calculated  tip  velocity  V and  tip  radius  R for  different  mesh  spacings  and 
time  steps.  The  results  were  computed  for  parameters  e = 0.005,  7 = 0.01,  A = 0.5, 
a = 400,  and  m — 0.05  on  a domain  with  Xl  — 2.0  and  Yl  — 1.0. 

e equal  to  0.005  or  smaller  the  interfacial  region  remains  thin.  We  conducted  a series  of 
calculations  in  which  only  the  interface  thickness  was  varied.  In  Table  3,  we  show  the 
dependence  of  the  dimensionless  tip  velocity,  radius,  and  temperature  on  the  interfacial 
thickness.  We  also  list  the  percentage  difference  between  the  computed  tip  temperature 
and  that  predicted  by  the  modified  Gibbs-Thomson  condition  Eq.  (16)  using  the  computed 
tip  velocity  and  radius;  this  provides  a quantitative  measure  of  the  accuracy  of  the  phase 
field  approximation  to  the  sharp  interface  model  Eq.  (14),  Eq.  (15)  and  Eq.  (24).  It  is 
evident  that  the  operating  state  of  the  dendrite  is  sensitive  to  the  interface  thickness.  How- 
ever, as  the  interface  width  is  reduced  the  error  in  the  modified  Gibbs-Thomson  equation  is 
diminished.  The  interface  moves  more  slowly  at  thinner  interface  widths.  Similar  behavior 
was  observed  in  the  spherically  symmetric  calculations  discussed  Section  3.3.  It  is  clear 
that  in  order  to  employ  the  phase  field  method  as  an  accurate  computational  approach  for 
approximating  the  solution  to  the  sharp  interface  model,  sufficiently  thin  interfaces  must 
be  taken  which  must  be  adequately  resolved  by  the  computational  mesh.  This  we  believe 
is  the  major  computational  issue  to  be  addressed  in  the  subsequent  development  of  the 
phase  field  method  as  an  accurate  computational  technique. 

From  the  above  discussion  in  order  to  adequately  resolve  the  interfacial  layers  we 
require  XX  — e,  which  for  a value  of  e of  2.5  xl0“^  produces  the  smallest  error  obtained 
in  the  phase  field  approximation.  However,  for  the  finite  difference  algorithm  with  a 
uniform  spatial  mesh  this  requires  an  impractical  amount  of  computing  resources  for  the 
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e 

V 

U{W^) 

X 100 

5.0  xlO-^ 

5.43 

6.79 

-3.86 

82  % 

3.3  xlO-^ 

3.89 

7.94 

-4.32 

19  % 

2.5  xlO-^ 

3.32 

10.23 

-3.90 

11  % 

Table  3:  The  calculated  tip  velocity  V,  tip  radius  and  interface  temperature  U for 
three  values  of  the  interface  thickness.  The  last  column  is  the  percentage  difference  of  the 
calculated  tip  temperature  to  the  temperature  Ugt  computed  from  the  modified  Gibbs- 
Thomson  Eq.  (24)  using  the  V and  R values.  The  results  were  computed  for  parameters 
7 = 0.01,  A = 0.5,  a = 400,  and  m = 0.05  on  a domain  with  Xl  = 2.0  and  = 1.0;  in 
each  case  the  mesh  spacing  was  given  by  AX  = e. 


simulation  of  dendrite  growth.  The  differences  given  in  Table  3 represent  a worst  case  and 
we  found  that  with  e = 5.0  xl0“^  the  error  in  satisfying  the  modified  Gibbs-Thomson 
equation  was  typically  25%.  We  go  on  to  describe  results  of  dendritic  growth  which  must 
be  regarded  in  this  context.  Our  aim  is  to  further  investigate  the  phase  field  model, 
and,  in  particular,  the  relation  of  the  computed  dendritic  solutions  to  existing  theories  of 
dendrites,  with  which  we  find  surprisingly  good  agreement. 

The  dimensions  of  the  computational  domain,  Xl  and  Yl  were  4.5  and  2.25,  respec- 
tively. Initially  the  solid  region  was  an  elHpse  with  semi-major  axis  Xo  = 0.5  and  semi- 
minor axis  2/o  ==  0.05;  both  the  solid  and  hquid  had  a dimensionless  temperature  of  —1.  We 
employed  a mesh  of  900  and  450  uniform  intervals  in  the  x and  y directions  respectively, 
corresponding  to  equal  mesh  sizes  AX  = AT  0.02,  with  a time  step  At  = lO""^.  Simple 
Euler  time  stepping  was  used  to  advance  the  solution  in  time.  This  mesh  was  the  finest 
we  could  practically  employ  consistent  with  doing  sufficient  runs  to  adequately  cover  the 
(7,  A)  parameter  space. 

Our  first  comparison  is  against  the  Ivantsov  similarity  solution  for  a parabohc  interface 
propagating  with  constant  velocity  in  the  direction  of  its  axis  of  symmetry  into  an  infinite 
undercooled  melt.  This  solution  assumes  that  the  interface  temperature  is  constant  and 
equal  to  the  melting  point,  and  so  interface  surface  energy  and  interface  kinetic  effects 
are  neglected.  It  predicts  in  two-dimensions  that  the  Peclet  number  is  related  to  the 
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dimensionless  undercooling  A as 


A = Vt?P  exp('P)  (25) 

where  V — 2vr  j k.  and  A = c[Tm  ~ ^oo)/-^-  Here  v and  r are  the  dimensional  tip  veloc- 
ity and  radius,  Tm  is  the  melting  temperature  which  is  equal  to  the  isothermal  interface 
temperature,  and  Too  is  the  far-held  temperature.  Our  phase  held  model  includes  both 
surface  energy  and  interface  kinetic  effects  and  so  the  interface  is  consequently  not  isother- 
mal with  a temperature  equal  to  the  melting  temperature.  This  is  conhrmed  in  Figure  4 
which  shows  that  the  dendrite  tip  temperature  U is  depressed  beneath  the  melting  point 
(given  as  zero  in  our  non-dimensionalization).  For  purposes  of  comparison  to  the  Ivantsov 
model,  which  only  accounts  for  conservation  of  heat,  we  estimate  the  undercooling  pa- 
rameter by  A = c{Ttip  — Ttnit)ILj  (=  U + 1),  where  Timt  is  the  initial  temperature  of 
the  hquid  in  the  phase  held  calculations.  The  initial  temperature  Tmit  is  a good  approx- 
imation to  Too  in  the  Ivantsov  model,  because  at  the  large  undercooHngs  we  employ  the 
temperature  gradient  is  weU-conhned  to  the  vicinity  of  the  dendrite  during  the  period  of 
constant  velocity  tip  growth.  In  Figure  5 we  compare  the  Peclet  numbers  estimated  from 
our  computations,  V = 2RV,  for  different  values  of  the  dimensionless  undercooling.  A, 
against  those  predicted  by  the  Ivantsov  formula,  Vj{A),  shown  as  a sohd  curve.  Here, 
the  function,  Vi{A),  is  obtained  by  inverting  the  Ivantsov  formula  Eq.  (25).  The  value 
of  the  anisotropy  parameter  is  7 — 10“^.  The  agreement  between  the  two  is  fairly  good, 
improving  as  the  undercooling  decreases.  In  Figure  6 we  plot  the  ratio  of  the  estimated 
Peclet  number  to  that  given  by  the  Ivantsov  formula  against  the  anisotropy  parameter 
for  A = 0.4.  We  chose  to  plot  the  ratio  V/Vi{A)j  because  the  estimated  undercoohng 
parameter,  and  hence  the  Peclet  number,  predicted  from  the  Ivantsov  theory  T/( A)  varies 
between  the  calculations  for  different  anisotropy  parameters.  We  see  from  Figure  6 that 
the  agreement  to  the  Ivantsov  theory  improves  as  the  anisotropy  diminishes.  It  would 
appear  that  extrapolation  of  the  almost  linear  dependence  of  V jVi^A)  to  7 = 0 would 
reveal  that  the  Ivantsov  solution  would  not  be  recovered  for  zero  anisotropy.  However,  as 
noted  above  no  steady  solution  was  found  for  7 < 10“^,  and  so  it  is  not  clear  whether 
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Undercooling  parameter,  A 


Figure  5:  The  Peclet  number  from  our  calculations  V against  the  estimated  undercooling 
parameter  A for  nickel.  The  sohd  circles  represent  runs  on  the  domain  Xl  — 4.5  and 
Yl  = 2.25  with  mesh  size  5 x 10“^;  the  open  squares  on  the  domain  Xl  = 2.0  and 
Yl  = 1.0  with  mesh  size  2.5  x 10“^.  The  solid  curve  is  the  prediction  from  the  Ivantsov 
formula  Eq.  (25). 
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Figure  6:  The  ratio  of  the  estimated  Peclet  number  V to  that  predicted  by  the  Ivantsov 
solution  Vi{K)  plotted  against  the  anisotropy  parameter  7 with  undercooling  parameter 


A = 0.4 
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there  is  any  justification  for  such  an  extrapolation.  In  Figure  7 we  compare  the  results  of 
our  computations  to  marginal  stability  theory.  Marginal  stabihty  theory  predicts  that 

O’"  - = 0.0192.  (26) 

Lvr^ 

In  this  figure  we  plot  the  values  of  a*  from  our  computations  against  the  dimensionless 
undercooling  parameter  A for  the  smallest  value  of  the  anisotropy  parameter  (7  = 10“^) 
for  which  we  obtained  a well  defined  steady  dendrite  tip  operating  state.  The  sohd  hori- 
zontal line  indicates  the  marginal  stability  result  given  above.  The  results  show  increasing 
agreement  with  the  marginal  stability  result  with  decreasing  undercooling.  The  depen- 
dence of  a*  on  A is  unexpected,  and  may  be  due  to  the  inclusion  of  interface  kinetics  in  the 
phase-field  model,  which  we  expect  to  be  significant  at  the  large  values  of  the  dimensionless 
undercooling  used  in  the  calculations,  but  which  are  absent  in  both  marginal  stability  and 
microscopic  solvabihty  theory.  In  Figure  8 we  display  against  the  anisotropy  parameter 
7 for  a fixed  value  of  A.  The  dashed  line  indicates  the  best  power  law  fit  through  the  data 
which  corresponds  to  a'*'  oc  7^-®^®.  Microscopic  solvability  theory  predicts  the  exponent  in 
the  power  law  to  be  1.75.  The  data  points  for  values  of  7 in  excess  of  2 x 10“^  consist 
of  dendrites  whose  tip  radius  [R  < 0.02)  is  well  less  than  the  interface  thickness  0.2  ) 
and  therefore  we  may  expect  these  results  to  be  less  reliable.  Also  shown  in  Figure  8 
as  solid  squares  are  the  results  of  microscopic  solvability  theory  given  by  Ben  Amar  [19] 
who  numerically  computed  the  value  of  cr*.  The  level  of  agreement  is  surprisingly  good, 
particularly  in  view  of  the  fact  the  we  are  conducting  computations  on  a finite  domain 
using  the  phase  field  model  with  a relatively  large  value  of  the  interface  width  which,  as 
we  have  discussed  above,  is  not  a particularly  good  approximation  to  the  classical  free 
boundary  problem  obtained  as  e — ^ 0.  Moreover,  this  free  boundary  problem  includes 
interface  kinetics  which  is  absent  in  the  microscopic  solvability  theory.  In  view  of  these 
limitations,  the  estimated  exponent  provides  some  evidence  to  suggest  that  microscopic 
solvabihty  theory  and  the  phase  field  model  in  two  dimensions  contain  a measure  of  agree- 
ment. The  experiments  conducted  by  Willnecker  et  al.  [20]  on  the  solidification  of  nickel 
measured  the  dendrite  tip  velocity  to  be  in  the  range  37  to  48  ms~^,  our  calculations 
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Figure  7;  The  calculated  value  of  a*  plotted  against  the  undercooHng  parameter  A.  The 
sohd  circles  represent  calculations  on  the  domain  Xl  = 4.5  and  Yl  = 2.25  with  mesh 
size  5 X 10~^;  the  open  squares  where  computed  on  the  domain  Xl  = 2.0  and  Yl  = l.O 
with  mesh  size  2.5  x 10“^.  The  horizontal  hne  corresponds  to  the  value  a*  = 0.0192  from 
marginal  stability  theory. 
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Figure  8:  The  calculated  value  of  a*  plotted  against  the  anisotropy  parameter,  7,  with 
undercooling  parameter  A = 0.4.  The  results  were  computed  on  the  domain  Xl  = 4.5 
and  Yl  = 2.25  with  mesh  size  5 x 10“^;  the  open  circles  represent  calculations  for  which 
a steady  operating  state  was  not  obtained.  The  sohd  squares  correspond  to  values  of  a* 
computed  by  Ben  Amar  [19]  on  the  basis  of  microscopic  solvability  theory.  The  dashed 
line  is  the  power  law  fit  to  the  solid  circles. 
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predict  tip  velocities  in  the  range  26  to  50ms for  the  anisotropy  parameter  in  the  range 
0.011  to  0.04.  The  value  of  7 for  most  materials,  including  nickel,  is  difficult  to  measure 
and  unknown.  However,  despite  this  uncertainty  in  the  value  of  7 our  calculations  reveal 
growth  rates  in  the  correct  range,  although  we  note  that  the  experiments  measured  three 
dimensional  dendrites  compared  to  our  two  dimensional  calculations. 

Finally,  we  use  the  present  phase  field  model  to  simulate  the  development  of  a side- 
branched  dendritic  structure.  We  found  that  the  formation  of  side-branches  was  depen- 
dent on  the  spatial  resolution  used.  For  spatially  underresolved  calculations,  with  ap- 
proximately 3-5  mesh  points  within  the  interfacial  layer,  uniformly  spaced  side-branches 
evolved.  However,  when  the  spatial  resolution  was  increased  the  side-branching  disap- 
peared. We  attribute  the  formation  of  side-branched  structure  on  coarser  meshes  to  the 
noise  associated  with  the  larger  truncation  error.  To  investigate  this  further,  we  introduced 
random  noise  into  the  phase  field  equation  for  computations  on  finer  meshes.  In  particular, 
we  added  the  term  AnTn  into  the  expression  in  square  brackets  in  the  discretized  form  of 
Eq.  (23),  where  An  is  the  amplitude  of  the  noise  and  is  a random  number  in  the  interval 
[—0.5, 0.5].  Physically,  this  represents  thermal  noise  at  the  interface  and  was  also  used  by 
Kobayashi  [1].  The  production  of  side  arms  was  very  sensitive  to  the  noise;  side  arms 
were  stimulated  for  values  of  An  as  low  as  2.5  xl0“^.  In  contrast  the  operating  state  of 
the  tip  was  insensitive  to  noise,  the  time  average  of  the  tip  velocity  and  radius  being  only 
very  weakly  affected.  This  provides  some  evidence  to  suggest  that  side  arm  formation  is 
a process  distinct  from  the  dynamics  of  the  dendrite  tip. 

In  Figure  9 we  show  a computation  that  displays  side-branch  formation.  The  phase 
field  and  isotherms  are  shown  at  three  different  dimensionless  times.  In  the  phase  field 
plot,  the  width  of  the  curve  represents  the  interval  0.1  < (j)  < 0.9,  and  thus  gives  an 
indication  of  the  thickness  of  the  interface.  For  this  computation,  we  again  employed  the 
parameter  values  corresponding  to  nickel  used  throughout  this  study.  For  the  simulation 
the  domain  has  dimensions  Xl  = 6.75  and  Yl  = 3.375.  The  mesh  spacing  was  AX  = 
AY  = 0.075,  and  the  time  step  was  At  = 5.0  xl0“^.  The  computation  starts  from  an 
elliptical  solid  region  with  semi-major  axis  Xq  = 0.15  and  semi-minor  axis  7/0  = 0.075.  The 
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dimensionless  parameters  are  e = 0.005,  7 = 0.01,  and  A = 0.5. 

5.  Conclusions 

We  have  conducted  an  extensive  detailed  study  of  the  phase  field  method  as  a computa- 
tional procedure  for  the  solution  of  free  boundary  problems  associated  with  sohdification, 
and  in  particular  dendritic  interfaces.  On  the  basis  of  the  simple  finite  difference  methods 
employed  here  we  find  that  in  simple  one  dimensional  geometries  accurate  solutions  can 
be  obtained.  However,  in  more  reahstic  two  dimensional  geometries  accurate  computa- 
tions require  interfaces  so  thin  that  the  need  to  accurately  resolve  them  require  computing 
resources  at  the  limit  of  current  computer  technology.  However,  less  accurate  computa- 
tions on  thicker  interfaces  provides  surprisingly  good  quantitative  agreement  with  both 
the  Ivantsov  and  microscopic  solvability  theories  of  dendrite  tip  growth.  Side  branch  for- 
mation in  these  computations  is  caused  by  the  introduction  of  very  small  amplitude  noise, 
which  leaves  the  operating  state  of  the  dendrite  tip  largely  unaffected. 

Because  of  the  need  to  use  very  fine  interface  thicknesses  we  beheve  that  further  use 
of  the  phase  field  model  as  an  accurate  computational  tool  for  the  computation  of  two 
and  three  dimensional  sohd  shapes  wiU  require  more  sophisticated  numerical  algorithms, 
possibly  employing  adaptive  finite  element  techniques. 
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Figure  9:  A simulation  of  a nickel  dendrite  with  the  phase  field  (on  the  left)  and  the 
temperature  field  (on  the  right)  shown  at  three  different  dimensionless  times  (r).  The 
dimensionless  parameters  are  e = 0.005,  7 = 0.01,  and  A = 0.5.  Random  noise  with  a 
one-percent  amphtude  was  introduced  into  the  calculation  to  stimulate  the  development  of 
the  side- branched  structure.  The  computational  domain  was  reflected  about  the  horizontal 
centerhne. 
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